package acewem.utilities.optimisation;

//import amesmarket.Support;
import acewem.market.ACEWEMmodel;
import acewem.market.GenCo;
import acewem.market.LSE;
import cern.colt.matrix.*;
import cern.colt.matrix.impl.*;
import java.util.*;

public class OptimalPowerFlow {
	
// Indices for Node Data parameters: form of {number of buses, Penalty coefficient}												 
	  private static final int NInd =       0;  				
	  private static final int PenCoefInd = 1;
	  
	  
// Indices for Branch Data parameters: form of {From, To, Branch Capacity, Reactance}
	  private static final int FromInd =  0;
	  private static final int ToInd =    1;
	  private static final int BCapInd =  2;
	  private static final int ReactInd = 3;
	  
	  
// Indices for Supply Offer parameters: form of {a,b,CapL,CapU}
	  private static final int aInd     = 0;
	  private static final int bInd     = 1;
	  private static final int CapLInd  = 2;
	  private static final int CapUInd  = 3;
	  														/** Number of Nodes */
	  private int K;
	  
	  private double KK;
	  														/** Number of Branches */
	  private int N;              
	  														/** Number of GenCos */
	  private int I;           
	  														/** Number of LSEs */
	  private int J;   
	  														/** GenCo's Data: {Gen ID, at Node, FCost, aT, bT, CapTL, CapTU, Initial Money in pounds} */
	 // private double[][] GenCoData; 
	  														/** LSE's Data: {LSE ID, at Node, L(0), L(1),L(2),L(3),L(4), ...,L(23)} */
	 // private double[][] LSEData; 	  
	  														/**  Node Data Array: {number of buses, Penalty coefficient} */
	  private double[][] NodeDataA;        
	  														/** Branch Data Array: {From, To, Branch Capacity, Reactance} */	  
	  private double[][] BranchDataA;
	  														/** Branch Data Array in SI units: {From, To, Branch Capacity, Reactance} */
	  private double[][] BranchDataASI;
	  														/** Hourly loadProfile Array = fixed demand : { each LSE demand over current hour } */	  
	  private double[]   LoadProfileA; 
	  														/** Hourly loadProfile Array in SI units = fixed demand : { each LSE demand over current hour } */
	  private double[]   LoadProfileASI;
	  														/** True Supply Offer Array: {a,b,CapL,CapU} */  
	  private double[][] SupplyOfferA; 
	  														/** True Supply Offer Array in SI units: {a,b,CapL,CapU} */ 
	  private double[][] SupplyOfferASI;
	  														/** Array of Generators' locations {gen1 at node 1, gen2 at node3,.., geni at nodej} */
	  private int[] GenAtNodeA;            
	  														/** Array of LSE locations {LSE1 at node 2, LSE2 at node5,.., LSEi at nodej} */
	  private int[] LSEAtNodeA;            	  												                 	  			                 
	  														/** Matrix to hold Branch Data:  {From, To, Branch Capacity, Reactance} */
	  private DoubleMatrix2D BranchDataM; 
	  														/** Matrix to hold Node Data: {number of buses, Penalty coefficient}  */
	  private DoubleMatrix2D NodeDataM;  
	  
	  
// Input for QuadProgJ (G,a,Ceq,beq,Ciq,biq)
	  														/** Matrix G =blockDiag [U  Wrr] of (I + K − 1)× (I + K − 1) dimension <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private DoubleMatrix2D GM;
	  														/** Vector a=[A1 ... Ai 0 ... 0] of (I + K − 1) dim, where Ai are a's from trues supply offer/pu <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private DoubleMatrix1D aM;
	  														/** The equality constraint matrix Ceq =[ II −Br ] of K × (I + K − 1) dim  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20)*/ 
	  private DoubleMatrix2D CeqM;
	  														/** The associated equality constraint vector  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */ 
	  private DoubleMatrix1D beqM;
	  														/** The inequality constraint constraint matrix <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20)*/ 
	  private DoubleMatrix2D CiqM;
	  														/** The associated inequality constraint vector <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 21)*/ 
	  private DoubleMatrix1D biqM;
	 
	  private QuadProgJ qpj;
	  
// Solution from QuadProgJ
	  														/** Power production quantity in MWs */
	  private double[] Commitment; 
	  														/** Voltage angle in radians */
	  private double[] VoltAngleR;  
	  														/** Voltage angle in degrees */
	  private double[] VoltAngleD;
	  														/** Full Voltage Angle including delta_1 = 0*/
	  double[] FullVoltAngle;
	  														/** LMP (locational marginal prices) == eqMultiplier */
	  private double[] LMP;  
	  														/** Lagrange multiplier corresponding to the inequality constraints */
	  private double[] IneqMultiplier;
	  														/** MinTVC in $/h = SUM (Ai*PGi + Bi*PGi^2) */
	  private double MinTVC; 
	  														/**  branchFlow_km in MWs = (1/x_km)(delta_k - delta_m) */
	  private double[] BranchFlow; 
	  														/** SUM(delta_k - delta_m)^2 */
	  private double sumSquaredAngleDifference;  
	  
	  
// Intermediate input to form (G,a,Ceq,beq,Ciq,biq)
	  														/** Vector B=[B1 ... Bi] of Ix1 dim, where Bi are b's from trues supply offer/pu */
	  private DoubleMatrix1D B;  
	  														/** An attribute matrix U = diag[2B1 , 2B2 , ··· , 2Bi] of Ix1 dim <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17)*/
	  private DoubleMatrix2D UM;    
	  														/** Reduced weight matrix  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17)*/
	  private DoubleMatrix2D WrrM; 
	  														/** Vector A=[A1 ... Ai] of Ix1 dim, where Ai are a's from trues supply offer/pu <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
	  private DoubleMatrix1D A;  
	  														/** II Matrix represents K Nodes(rows)  and I Generators(columns) located on them */
	  private DoubleMatrix2D IIM;
	   														/** Nx1 matrix  of branches' capacities  */
	  private DoubleMatrix1D BCapM;   
	  														/** Ix1  matrix of generator's true lower capacity */
	  private DoubleMatrix1D CapTLM;  
	  														/** Ix1 matrix of generator's true upper capacity */
	  private DoubleMatrix1D CapTUM; 
	  														/** (N×I) zero O matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */ 
	  private DoubleMatrix2D OtM;
	  														/** (I x K-1) zero O matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21)  */
	  private DoubleMatrix2D OpM;
	  														/** I×I Identity Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	  private DoubleMatrix2D IpM;
	  														/** (NxN)Diagonal Admittance Matrix, where Bkm=1/Xkm <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	  private DoubleMatrix2D DAM;
	  														/** (N x K-1) Reduced Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19)  */
	  private DoubleMatrix2D RAM;
															/** True Supply Offer Matrix: {a,b,CapL,CapU} */
	  private DoubleMatrix2D SupplyOfferMatrix; 														
	  														/** (KxK) Forms Voltage Angle Difference Weight Matrix <br>@see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */
	  private double[][] VADWM;               
	  														/** (K-1)x(K-1) Reduced Voltage Angle Difference Weight Matrix <br>@see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17) */
	  private double[][] RVADWM;             
	  														/** (Nx2) Matrix which holds the data for Branches' connections, e.g.  {{1,2},{1,4},{2,3},...} */
	  private double[][] BranchIndex;            
	  														/** (KxK) Negative Susceptance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
	  private double[][] NSM;     
	  														/** (KxK) Bus Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private double[][] BAM;       
	  														/** Reduced Bus Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */
	  private double[][] RBAM;
	  														/**(NxK) Adjacency Matrix <br>@see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19)  */
	  private double[][] AM;         
	  														/** Temp matrix to form DA matrix*/
	  private double[][] dam;  
	  														/** Temp matrix to form RA matrix */
	  private double[][] ram;       
	  														/** A factory producing dense matrices 1D */
	 
	  private DoubleFactory1D Matrix1d = DoubleFactory1D.dense; 
	  														/**A factory producing dense matrices 2D*/
	  private DoubleFactory2D Matrix2d = DoubleFactory2D.dense; 
	  														/** Func is for using Colt's methods e.g. diagonal(), identity(), etc. */
	  private cern.jet.math.Functions Func = cern.jet.math.Functions.functions;

	 
	  private ACEWEMmodel market;
	  
	  Hashtable<String, Double> lmp = new Hashtable<String, Double>();
	  
//C-O-N-S-T-R-U-C-T-O-R------------------------------------------------------------------------------------------------------------------------------  
	  public OptimalPowerFlow(ACEWEMmodel mkt)
	  {
		  this.market=mkt;
		//  this.hour=h;

        	I      = mkt.genCoList.size();
        	J      = mkt.lseList.size();
        	KK      = mkt.nodeData.get("buses");
        	K = (int)KK;
        	N 	   = mkt.branchList.size();
        	
        	
        	//System.out.println(I);
        	//System.out.println(J);
        	//System.out.println(K);
        	//System.out.println(N);

        	
        	//init            = new Initials.InitParameters();
        	NodeDataA    	= new double[1][2];       	
        	BranchDataA	    = new double[N][4];
        	BranchDataASI   = new double[N][4];        	
        	LoadProfileA 	= new double[J];  
        	LoadProfileASI	= new double[J];         	
        	SupplyOfferA 	= new double[I][4];
        	SupplyOfferASI  = new double[I][4];
        	BranchIndex     = new double[N][2];
        	//GenCoData       = new double[I][8];
        	//LSEData         = new double[J][26];
        	GenAtNodeA 		= new int[I];	
        	LSEAtNodeA 		= new int[J];
        	

	    	

		}
	  
//------------------------------------------------------------------------------------------------------------------------------------------    	

	  public void solveHourlyPowerFlows()
	  {
		  for (int hour=0; hour<market.Hours; hour++)
		  {	
			//System.out.println("solveHourlyPowerFlows"+ hour); 
			buildDataForOPF(hour);
	    	convertSItoPU();      
	    	solveOPF(hour);
		  }
		  
	  } 
	  
	  
//------------------------------------------------------------------------------------------------------------------------------------------    	
	  private void buildDataForOPF(int hour){
		  
		//Nodes----------------------------------------------------------------------------       	       	
      	NodeDataA[0][0]    = market.nodeData.get("buses");
      	NodeDataA[0][1]    = market.nodeData.get("penalty");
      	                	       
//Branches-------------------------------------------------------------------------       	
      	for (int n=0; n<N; n++)
      	{
      		int id = n+1;  		
      		Hashtable<String, Double> BranchHash = new Hashtable<String, Double>();
      		BranchHash =(Hashtable) market.branchList.get("branch"+id);
      		    		
      		BranchDataA[n][0] = ((Double)BranchHash.get("from")).doubleValue();
      		BranchDataA[n][1] = ((Double)BranchHash.get("to")).doubleValue();
      		BranchDataA[n][2] = ((Double)BranchHash.get("capacity")).doubleValue();
      		BranchDataA[n][3] = ((Double)BranchHash.get("reactance")).doubleValue();       	        		        		
      	} 
      	
//LSEdemand-------------------------------------------------------------------------        
      	for (int j=0; j<J; j++)
      	{
      		int id = j+1;
      		LoadProfileA[j] = market.lseList.get("lse"+id).electricityDemand.get(Integer.toString(hour)); 
      		
				//if (id ==2){
				//		System.out.println("222opfXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx");
				//		System.out.println("hour= "+hour+ "  "+LoadProfileA[j]);
				//}
      	}
      //	LoadProfileA = Support.correctRoundingError(LoadProfileA);
//GenCo------------------------------------------------------------------------------       
      	for (int i=0; i<I; i++)
      	{    		
      		int id = i+1;
      		GenCo gen = market.genCoList.get("genco"+id);
      		
      	/*	SupplyOfferA[i][0] = gen.reportedSupplyOffer.get("aR");
      		SupplyOfferA[i][1] = gen.reportedSupplyOffer.get("bR");
      		SupplyOfferA[i][2] = gen.reportedSupplyOffer.get("capRL");
      		SupplyOfferA[i][3] = gen.reportedSupplyOffer.get("capRU"); 
      		*/
      		SupplyOfferA[i][0] = gen.getReportedSupplyOffer().get("aR");
      		SupplyOfferA[i][1] = gen.getReportedSupplyOffer().get("bR");
      		SupplyOfferA[i][2] = gen.getReportedSupplyOffer().get("capRL");
      		SupplyOfferA[i][3] = gen.getReportedSupplyOffer().get("capRU"); 
      		
      		
      		/*	
      		    SupplyOfferA[i][0] = gen.supplyOffer.get("aT");;
          		SupplyOfferA[i][1] = gen.supplyOffer.get("bT");
          		SupplyOfferA[i][2] = gen.supplyOffer.get("capTL");
          		SupplyOfferA[i][3] = gen.supplyOffer.get("capTU"); 
      		*/
      	//	System.out.println("id="+id+"  "+ SupplyOfferA[i][0]+"  "+SupplyOfferA[i][1]+"   "+SupplyOfferA[i][2]+"  "+SupplyOfferA[i][3]);
      	}

      	//SupplyOfferA  = supplyoffer;
      	//GenCoData     = gencodata;
      	//LSEData       = lsedata;       	

      	
      	GenCo gen=null;
	    	for (int g=0; g<I; g++ ){
	    		int id = g+1;
	    		gen = market.genCoList.get("genco"+id);
	    		GenAtNodeA[g] = gen.getNode();
	    	}
	    
	    	LSE lse=null;
	    	for (int l=0; l<J; l++ ){
	    		//LSEAtNodeA[l] = (int)LSEData[l][1];
	    		int id = l+1;
	    		lse  = market.lseList.get("lse"+id);	    		
	    		LSEAtNodeA[l] = lse.node;
	    	}		  
	  }
	  
//------------------------------------------------------------------------------------------------------------------------------------------
	/** Converts Branch Capacity,Reactance, LSE load profiles, */   
	 private void convertSItoPU(){
		 
			   acewem.initials.InitParameters init  =   new  acewem.initials.InitParameters();
			   
		   for(int n=0; n<BranchDataA.length; n++){
		        
		    	BranchDataASI[n][0] = BranchDataA[n][0];
		    	BranchDataASI[n][1] = BranchDataA[n][1];
		    	
	// Convert Branch Capacity from SI to PU
		    	BranchDataASI[n][2] = BranchDataA[n][2]/init.getBaseS();
		    	
	// Convert Reactance from SI to PU, x(pu) = x/Zo = x/(Vo^2/So) = (x*So)/Vo^2
		    	BranchDataASI[n][3] = (BranchDataA[n][3]*init.getBaseS())/(init.getBaseV()*init.getBaseV());
		      }
		        //BranchDataA = Support.correctRoundingError(BranchDataA);
			
	// SI to PU conversion for supply offer and load profile
		      double sumCapL=0.0;
		      double sumCapU=0.0;
		      
	// Convert hourly Load Profile from SI to PU
		      for(int n=0; n<J; n++){
		    	  LoadProfileASI[n] = LoadProfileA[n]/init.getBaseS();}
	//  LoadProfileA = Support.correctRoundingError(LoadProfileA);
		      
		      
		     for(int i=0; i<I; i++){
		         sumCapL+=SupplyOfferA[i][CapLInd];
		         sumCapU+=SupplyOfferA[i][CapUInd];
		         //System.out.println("\n supplyOfferByGen-> dMaxGenCapacity "+supplyOfferByGen[i][CAP_MAX]+" for GenCo :"+i);

	// Convert A from SI to PU-adjusted
		         SupplyOfferASI[i][aInd] = SupplyOfferA[i][aInd]*init.getBaseS();

	// Convert B from SI to PU-adjusted
		         SupplyOfferASI[i][bInd] = SupplyOfferA[i][bInd]*init.getBaseS()*init.getBaseS();

	// Convert CapMin from SI to PU
		         SupplyOfferASI[i][CapLInd] = SupplyOfferA[i][CapLInd]/init.getBaseS();

	// Convert CapMax from SI to PU
		         SupplyOfferASI[i][CapUInd] = SupplyOfferA[i][CapUInd]/init.getBaseS();
		         
		
		          
		     }
		     // TSupplyOfferArray = Support.correctRoundingError(TSupplyOfferArray);
		   //  SupplyOfferASI = Support.correctRoundingError(SupplyOfferASI);
			
			  SupplyOfferMatrix = new DenseDoubleMatrix2D(SupplyOfferASI);
			  
			//  System.out.println("SupplyOfferMatrix:" + SupplyOfferMatrix);
			  
			  
			  BranchIndex = new double[N][2];  //e.g., {{1,2},{1,4},{2,3},...}
			
		      NodeDataM = new DenseDoubleMatrix2D(NodeDataA);
		      BranchDataM = new DenseDoubleMatrix2D(BranchDataASI);
		      
		      //System.out.println("G: " + NodeDataM);
		   }	 
//------------------------------------------------------------------------------------------------------------------------------------------    	
    /** Calculates parameters for the Optimal Power Flow  */
    private void solveOPF(int hour){
		
		
	    setBranchIndex();
	    
	    formGM();
	    formaM();
	    formCeqM();
	    formbeqM();
	    formCiqM();
	    formbiqM();
		
/*	    System.out.println("G: " + GM);
	    System.out.println("a: " + aM);
	    System.out.println("Ceq': " + CeqM.viewDice());
	    System.out.println("beq: " + beqM);
	    System.out.println("Ciq': " + CiqM.viewDice());
	    System.out.println("biq: " + biqM); 
*/	
	    
	   // System.out.println("Ceq': " + CeqM.viewDice());
	    
	    qpj = new QuadProgJ(GM,aM,CeqM,beqM,CiqM,biqM);

	    boolean bHaveSolution=qpj.getIsFeasibleAndOptimal();

	    Commitment 	   = new double[I];   		    
	    VoltAngleR     = new double[K-1]; 			
	    VoltAngleD 	   = new double[K-1];           
	    IneqMultiplier = new double[2*N+2*I];
	    LMP            = new double[K];
	    MinTVC         = 0;  						
	    BranchFlow     = new double[N];          
	    FullVoltAngle  = new double[K];      
	    sumSquaredAngleDifference = 0;
	    
	    if(bHaveSolution) { 										// QuadProgJ has a solution
// OPF solution for (p_{G1},...,p_{GI}) in SI
	    	GenCo gen=null;
	        for(int i=0; i<I; i++){
	        	Commitment[i] = qpj.getMinX()[i]*market.getBaseS();
	        	int id = i+1;
	        	gen = market.genCoList.get("genco"+id);
	        	gen.getCommitment().put("commit"+hour, Commitment[i]);
	        	
					//if (i ==1){
						//System.out.println("305genXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx");
					//		System.out.println("hour= "+hour+ "  "+Commitment[i]);
				//	}
	        	
	        }	         
	      
// OPF solution for (delta_2,...,delta_K)
	        for(int k=I; k<I+K-1; k++){
	            VoltAngleR[k-I] = qpj.getMinX()[k];  
	         }   

// Convert voltage angle from radian to degree
	        for(int k=1; k<K-1; k++){
	            VoltAngleD[k] = (VoltAngleR[k]*180)/Math.PI;
	            //System.out.println("V Angle" + VoltAngleD[k]);	          
	        }

	        for(int j=0; j<2*N+2*I; j++){
		        IneqMultiplier[j] = qpj.getIneqMultipiers()[j]/market.getBaseS();
		        //System.out.println("IneqMultiplier" +  IneqMultiplier[j]);
     	    }
	        

	        for(int l=0; l<I; l++){
	        	MinTVC = MinTVC + (A.get(l)/market.getBaseS())*Commitment[l]
	            +(B.get(l)/(market.getBaseS()*market.getBaseS()))*Commitment[l]*Commitment[l];
	            //System.out.println("TVC" +  minTVC);
	        }

	        for(int k=1; k<K; k++){
	          FullVoltAngle[k] = VoltAngleR[k-1];
	        }
	        
	        for(int n=0; n<N; n++){
	            BranchFlow[n] = (1/BranchDataASI[n][3])*(FullVoltAngle[(int)BranchIndex[n][0]-1]
	                         - FullVoltAngle[(int)BranchIndex[n][1]-1])*market.getBaseS();

	            sumSquaredAngleDifference = sumSquaredAngleDifference
	              + Math.pow((FullVoltAngle[(int)BranchIndex[n][0]-1]
	                          - FullVoltAngle[(int)BranchIndex[n][1]-1]),2);
	            // System.out.println("BranchFlow" +  BranchFlow[n]);
	        }
	        
// LMP: locational marginal prices in SI
	        for(int k=0; k<K; k++){	        	       	        		        	
	          LMP[k] = qpj.getEqMultipliers()[k]/market.getBaseS();	        
	          int id = k+1;
	          lmp.put("lmp"+id, LMP[k]);
	          
			//	if (k ==1){
				//	System.out.println("360opfXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx");
				//		System.out.println("hour= "+hour+ "  "+LMP[k]);
			//	}
	          
	        }		       
	        market.lmpAtHour.put("lmpathour"+hour, lmp); 	        
        }
 
	  // System.err.println("Error: qpj.getIsFeasibleAndOptimal(): check the solveOPF() method");    
	}
	
//------------------------------------------------------------------------------------------------------------------------------------------    	
/** Forms Voltage Angle Difference Weight Matrix
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */	 	
	private void formGM() {
		
	    formVADWM();
	    formRVADWM();
	    
	     B = new DenseDoubleMatrix1D(SupplyOfferMatrix.viewColumn(bInd).toArray());         // Copies converted (/pu) values of b from supply offer into array B
	     
	     //System.out.println("B:" + B);
	     
	     UM = new DenseDoubleMatrix2D(Matrix2d.diagonal(B.assign(Func.mult(2))).toArray()); // Matrix U has 2*B diagonal and 0 for other elements
	     WrrM = new DenseDoubleMatrix2D(RVADWM);
	     GM = new DenseDoubleMatrix2D(I+K-1,I+K-1);
	     GM.assign(Matrix2d.composeDiagonal(UM,WrrM));
	 }
	 
	 
	 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Voltage Angle Difference Weight Matrix
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */ 
	 private void formVADWM(){
		 
		 VADWM = new double[K][K];
		 
		    for(int n=0; n<N; n++){
		      VADWM[(int)BranchIndex[n][0]-1][(int)BranchIndex[n][1]-1]
		          = -2*NodeDataA[0][PenCoefInd];  //NOTE: there should be a factor 2 in front of penaltyCoeff
		      VADWM[(int)BranchIndex[n][1]-1][(int)BranchIndex[n][0]-1]
		          = -2*NodeDataA[0][PenCoefInd];
		    }
		    for(int i=0; i<K; i++){
		      for(int j=0; j<K; j++){
		        if(j==i){
		          for(int k=0; k<K; k++){
		            if(k!=i){
		            	VADWM[i][j] = VADWM[i][j] - VADWM[i][k];
		            }
		          }
		        }
		      }
		    }
		 //   VADWMatrix = Support.correctRoundingError(VADWMatrix); // Need a support Class
		  }
	 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Reduced Voltage Angle Difference Weight Matrix by deleting of VADWMatrix first row and first column
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 17) */ 	
   private void formRVADWM(){
	   
	   RVADWM = new double[K-1][K-1];
	   
		    for(int i=0; i<K-1; i++){
		      for(int j=0; j<K-1; j++){
		    	  RVADWM[i][j] = VADWM[i+1][j+1];
		      }
		    }
		  }
   
   
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Matrix of a coefficients from the supply offer
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) 
*  alpha*P+beta*P^2
*  */ 
   private void formaM(){
     A = new DenseDoubleMatrix1D(SupplyOfferMatrix.viewColumn(aInd).toArray());

     aM = new DenseDoubleMatrix1D(I+K-1);
     aM.viewPart(0,I).assign(A); // viewPart(fromIndex, width)       // ?????????????????????????????
   }
   
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms  Adjacency Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */ 
      private void formAM(){  
    	  
    	  AM = new double[N][K];
    	  
    	    for(int n=0; n<N; n++){
    	      for(int k=0; k<K; k++){
    	        if(k == BranchIndex[n][0]-1){
    	        	AM[n][k] = 1;
    	        }
    	        else if (k == BranchIndex[n][1]-1){
    	        	AM[n][k] = -1;
    	        }
    	        else{
    	        	AM[n][k] = 0;
    	        }
    	      }
    	    }
    	  }
      
      
      
//------------------------------------------------------------------------------------------------------------------------------------------          
/** Forms Reduced  Adjacency Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */  
      private void formRAM(){
    	  
    	  ram = new double[N][K-1];
    	  
    	    DoubleMatrix2D amx = new DenseDoubleMatrix2D(AM);
    	    ram= amx.viewPart(0,1,N,K-1).toArray();
    	    RAM = new DenseDoubleMatrix2D(ram);
    	  }
      
 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms II Matrix which represents K Nodes(rows)  and I Generators(columns) located on them.
 * <f/> If generator Ii is located at node Ki then the matrix element has value of 1 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */    
   private void formIIM(){
	   
   IIM = new DenseDoubleMatrix2D(K,I);
   for(int k=0; k<K; k++){
     for(int i=0; i<I; i++){
       if(GenAtNodeA[i]==k+1){
         IIM.set(k,i,1);  // change set to setQuick() later
       }}}}

//------------------------------------------------------------------------------------------------------------------------------------------      
/** Forms Bus Admittance Matrix.
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */    
    private void formBAM(){
    	  
    	BAM = new double[K][K];
    	
    	  for(int i=0; i<K; i++){
    	      for(int j=0; j<K; j++){
    	        if(j==i){
    	          for(int k=0; k<K; k++){
    	            if(k!=i){
    	            	BAM[i][j] = BAM[i][j] + NSM[i][k];
    	            }
    	          }
    	        }
    	        else
    	        	BAM[i][j] = - NSM[i][j];
    	      }
    	    }
    	  
    	  
      }

//------------------------------------------------------------------------------------------------------------------------------------------       
/** Forms Reduced Bus Admittance Matrix.
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
    private void formRBAM(){
    	
    	RBAM = new double[K-1][K];
    	
    	    DoubleMatrix2D na = new DenseDoubleMatrix2D(BAM);
    	    RBAM = na.viewPart(1,0,K-1,K).toArray();
    	    }
   
//------------------------------------------------------------------------------------------------------------------------------------------        	    
 /** Forms Negative Susceptance Matrix.
  *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
    private void formNSM(){
    	
    	NSM = new double[K][K];
    	
    	for(int n=0; n<N; n++){
    		NSM[(int)BranchIndex[n][0]-1][(int)BranchIndex[n][1]-1]
    		          = 1/BranchDataM.viewColumn(ReactInd).toArray()[n];
    		NSM[(int)BranchIndex[n][1]-1][(int)BranchIndex[n][0]-1]
    		          = 1/BranchDataM.viewColumn(ReactInd).toArray()[n];
    		    }
    	//System.out.println("G: " + BranchDataM);
    	//for (int aa=0; aa<K; aa++){
    	//	for (int bb=0; bb<K; bb++){
    		//	System.out.print(1/NSM[aa][bb]+"  ");
    		//} 
    	////	System.out.println();
    	//}
     }
    	    
 //------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Equality Constraint Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */ 
   // CeqTranspose = (II, -Br'); Ceq = CeqTranspose'; where Br' is rBusAdm here
   private void formCeqM(){
	   
	    formNSM();
	    formBAM();
	    formRBAM();
	    formIIM();
       
	    DoubleMatrix2D rBAM = new DenseDoubleMatrix2D(RBAM);
        DoubleMatrix2D[][] parts = {{ IIM, rBAM.viewDice().assign(Func.neg)}};
        DoubleMatrix2D CeqTMatrix = new DenseDoubleMatrix2D(K,I+K-1);
        CeqTMatrix.assign(Matrix2d.compose(parts));
        CeqM = new DenseDoubleMatrix2D(I+K-1,K);
        CeqM.assign(CeqTMatrix.viewDice());
   
   }   
   
   
//------------------------------------------------------------------------------------------------------------------------------------------    

// FDemand = someFunction(atNodeByLSE, loadProfile); beq = FDemand
/** Forms an Equality Constraint Vector 
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */
       private void formbeqM(){

    	 beqM = new DenseDoubleMatrix1D(K);
         for(int k=0; k<K; k++){
           double L=0;
           for(int j=0; j<J; j++){
             if(LSEAtNodeA[j]==k+1){
               L = L + LoadProfileASI[j];
             }
           }
           beqM.set(k,L);
         }
     
   }   

//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms an Inequality Constraint Matrix Ciq 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */     

       private void formCiqM(){

    	   formOtM();
   	       formOpM();
   	       formDAM();
   	       formAM();
   	       formRAM();
   	       formIpM();
    	   
   	    
       
           DoubleMatrix2D[][] parts = {
           { OtM,    	   DAM.zMult(RAM,null)         }, 
           { OtM,             DAM.copy().assign(Func.neg).zMult(RAM,null)},  // (1) changed upside-down           
           { IpM,                             OpM                   },
           { IpM.copy().assign(Func.neg),     OpM                  }};
           
           
           DoubleMatrix2D CiqTMatrix = new DenseDoubleMatrix2D(2*N+2*I,I+K-1);
           CiqTMatrix.assign(Matrix2d.compose(parts));
           CiqM = new DenseDoubleMatrix2D(I+K-1,2*N+2*I);
           CiqM.assign(CiqTMatrix.viewDice()); 
       }
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms Diagonal Admittance Matrix, where Bkm=1/Xkm 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */       
       private void formDAM(){
    	   
    	   dam = new double[N][N];
    	   
    	   dam = Matrix2d.diagonal(BranchDataM.copy().viewColumn(ReactInd).assign(Func.inv)).toArray();
    	   DAM = new DenseDoubleMatrix2D(dam);   //NOTE: Have to keep .copy(), otherwise reactance will be 1/reactance
    	  
    	   
    	   
       	//System.out.println("G: " + BranchDataM);
       	for (int aa=0; aa<N; aa++){
       		for (int bb=0; bb<N; bb++){
       			System.out.print(dam[aa][bb]+"  ");
       		} 
       		System.out.println();
       	}
    	   
    	   
    	   
    	   
    	  
    	   }
    	   
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms O an N × I zero Matrix 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
       private void formOtM(){
    	   
    	   OtM   = new DenseDoubleMatrix2D(N,I);     	   
       }
 
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms O an (I x K-1) zero Matrix 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
        private void formOpM(){
           	   
        	OpM = new DenseDoubleMatrix2D(I,K-1); 
        }
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms I×I Identity Matrix 
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
        private void formIpM(){
                   	   
        	IpM = new DenseDoubleMatrix2D(I,I).assign(Matrix2d.identity(I));
        }  
        
        
 //------------------------------------------------------------------------------------------------------------------------------------------	 
	 private void setBranchIndex(){
		    BranchIndex = BranchDataM.viewPart(0,FromInd,N,2).toArray();
		    }
	 


//------------------------------------------------------------------------------------------------------------------------------------------
// biq = (-pU, -pU, capL, -capU)
/** Forms (2N+2I)×1 the associated inequality constraint vector 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	 private void formbiqM(){

	
	 
		 BCapM =new DenseDoubleMatrix1D(BranchDataM.viewColumn(BCapInd).toArray()); // !!! Check if it's working	 
		 CapTLM = new DenseDoubleMatrix1D(SupplyOfferMatrix.viewColumn(CapLInd).toArray());
		 CapTUM = new DenseDoubleMatrix1D(SupplyOfferMatrix.viewColumn(CapUInd).toArray());

	 
		 DoubleMatrix1D[] parts = {BCapM.copy().assign(Func.neg), BCapM.copy().assign(Func.neg),
    		                  CapTLM, CapTUM.copy().assign(Func.neg)};
    
    
		 biqM = new DenseDoubleMatrix1D(2*N+2*I);
		 biqM.assign(Matrix1d.make(parts));
    
     }
	 


}












